Green's function of a dressed particle 
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We present a new, highly efficient yet accurate approximation for the Green's functions of dressed 
particles, using the Holstein polaron as an example. Instead of summing a subclass of diagrams {e.g. 
the non-crossed ones, in the self-consistent Born approximation (SCBA)), we sum all the diagrams, 
but with each diagram averaged over its free propagators' momenta. The resulting Green's function 
satisfies exactly the first six spectral weight sum rules. All higher sum rules are satisfied with great 
accuracy, becoming asymptotically exact for coupling both much larger and much smaller than the 
free particle bandwidth. Possible generalizations to other models are also discussed. 

PACS numbers: 71.38.-k, 72.10.Di, 63.20.Kr 



One of the most fundamental problems in both high- 
energy and condensed matter physics is to understand 
what happens when a particle couples to an environment, 
in particular what are the properties of the resulting ob- 
ject, consisting of the bare particle dressed by a cloud of 
excitations. This type of problem arises again and again 
as couplings to new kinds of environments are studied. 

The most desirable quantity to know is the Green's 
function G{k,uj) of the dressed particle - its poles mark 
the cigenspectrum, while the associated residues contain 
information on the eigenfunctions. Moreover, the spec- 
tral weight A{k,LLj) — ~^lmG{k,uj) can be directly mea- 
sured experimentally using Angle-Resolved Photoemis- 
sion Spectroscopy jlJ. Recently, such work has reignited 
a debate on whether the carriers in high-Tc cuprates are 
polarons, that is, electrons dressed by phonons 

G{k, uj) is the sum of an infinite number of diagrams 
corresponding to an expansion to all orders in the cou- 
pling strength Diagrammatic Monte Carlo (DMC) 
can perform the numerical summation of all diagrams jj| . 
Other ways to find G{k,uj) are from exact diagonaliza- 
tions (ED) of small systems, variational methods, Den- 
sity Matrix Renormalization Group in one-dimension, 
etc IE IE S 111- However, these methods require con- 
siderable computational resources, are time consuming, 
and often limit themselves to finding only the low-energy 
properties, such as the ground-state energy. 

To our knowledge, there are only two easy-to-estimate 
approximations for G{k,uj). One is the the SCBA which 
consists in summing only the non-crossed diagrams. Be- 
cause the percentage of diagrams kept decreases fast 
with increasing order, SCBA fails badly at strong cou- 
plings. The other approximation, obtained from a mod- 
ified Lang-Firsov (MLF) approach |^, is exact both for 



E(fc,w) 



FIG. 1: Diagrammatics for E, where G = Go + GqSG. 




zero coupling and for zero bandwidth, however results for 
finite bandwidth/coupling are rather poor (see below). 

In this Letter we find a new approximation for G{k, lu), 
which is as easy to estimate as SCBA and MLF, but is 
highly accurate over most of the parameter space. We val- 
idate this both by comparison against numerical results, 
and by investigating its sum rules. Most of the discussion 
here is limited to the Holstein model, for which many nu- 
merical results are available. Possible generalizations for 
other models are briefly discussed at the end. 

Consider the Holstein Hamiltonian: 

k k,q 

which contains the kinetic energy of the free particle 
(electron) , a branch of Einstein bosons (optical phonons) 
and the linear coupling between particle and bosons. 
When needed, we use e^: — — 2i ^^^^ cos(fcaa), for 
nearest-neighbor hopping on a d-dimensional cubic lat- 
tice with sites and lattice constant a. The spin of the 
particle and whether it is a boson or a fermion is irrel- 
evant. Sums over momenta are over the first Brillouin 
zone (BZ), -f < A:^, < f , a = l,d. 

The one-particle Green's function of a M-particle sys- 
tem is G(fc,T) = -i($M|T[cg(T)ct(0)]|$M), with |$Af) 

the M-particle ground-state, T the time ordering oper- 
ator, and cj(t) = exp(i7iT)cj: exp(— iTir) (we set h = 
1) l3|- The spectrum of a single dressed particle (po- 
laron) is obtained for M = 0. The ground state (GS) in 
the absence of particles is the vacuum |$o) = |0)i ^nd 
since 7i|0) — 0, the polaron's Green's function simplifies 
to: 



G{k,T) = -tQ{T){0\cf:, 
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where 8(t) is the Heavisidc function. From Eqs. and 
we derive the equation of motion: 



i—G{k, t) = 5{t) + e^G{k, r) + ^ V Fi (fc, q, , r) 

91 



2 



where Fi(fc,gi,r) = -ie(r)(0|cgexp(-iHr)ct_^ &t |0). 
Continuing in this vein, we generate an infinite hierarchy 
of coupied equations of motion for the Green's functions 
F„(fc,<fi, . . . ,g„,T) = -*e(T)(0|c^~e-*«-4__.^6j-.^ • • -fetjO), 

with qx = X]"=i ^ ^^"^ G{k,T) = Ft){k,T). In the fre- 
quency domain, these equations of motion become: 

G(fc, c.) = Go (fc, c.) [1 + ^ Fi (fc, gl , (3) 

91 

and for any ri > 1, 

Fn{k,qi, . . .,qn,uj) = -j=Go{k - qr.uj - uuje) 
V N 



'x[^Fn-i{k,qi, . . .,qi_i,qi+i, . . . ,g„,w) 
1=1 

+ ^ Fn+i{k,qi,...,qn,q„+i,uj)] (4) 

where Go{k,uj) = [tj — + i?;]^^ is the non-interacting 
one-particle Green's function. This system generates 
the expected diagrammatic expansion for the self-energy 
= Go ^(fc,w) - G-^{k,uj), shown in Fig. □ 
Let /„(fc, uj) = Fn iK gi, . . . , gk, c^). In 

terms of these, Eq. |j3Jl becomes G(fc, w) = Go(fc, u;)[l -I- 
gVNfi{k,U!)]. The equations for fn{k,Lu), n > 1, are 
obtained by summing Eqs. Q over all phonon momenta. 
Of the two terms on the right-hand side, the first one 
can be expressed in terms of fn~i{k,Lu) exactly, but the 
second one requires an approximation. We replace: 

'^Go{k ~ qT,uj ~ nuJE)Fn+i{k,qi, ■ ■ ■,qn+i,i^) 

^N^+^go{uJ~nLOE)fn+i{k.u:) (5) 



(ji,...,g„ + i 



where 



(6) 



The justification is that gV = Yl7=i 9j takes, with equal 
probability, any value in the first Brillouin zone. Replac- 
ing Go(/c - qT,uj~ uuje) {Go{k - qT,oj- nu}E))qT = 
gQ{uj~nujE) allows us to also write this term as a function 
of fn+i{k,u}) only. We discuss below the meaning of this 
momentum average (MA) in terms of diagrams; however, 
note that for hopping i = this MA approximation be- 
comes exact, because for i = all Green's functions are 
independent of momenta. This suggests that MA should 
be vaHd at least in the strong-couphng regime t/g <C 1. 
As we show later, its validity range is in fact much wider. 
With this approximation, Eqs. Q become /„(fc, a;) = 

go{LO - uuje) 



^f^_,{k,u;)+gVNU+i{k,Lu) 
cursive chain has a continued-fraction solution 



This re- 

The re- 
sulting Green's function can be cast in the usual form 



GMA{k,uj) ■ 



irj] ^ 



where: 



2g'^go{uj - ujE)goi^ - '^^e) 



(7) 
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As pointed already out, if i = 0, in which case goioj) = 
(aH-i?7)~^, this expression is exact. Indeed, one can show 
[loj that it equals the expected Lang-Firsov result Q 
[X={g/UE?]: 



Giio) 



E 

n=0 



X" 
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n\ UJ + Xlue — niUE + it] 



(8) 



This mapping was used before in a Dynamical Mean Field 
Theory (DMFT) study of this problem, which also pro- 
duces an approximation for the Green's function 0. In 
fact, Eq. Q looks similar to T,]j]\iFTi^)', however, our 
go{uj) is not a solution of the self-consistent DMFT equa- 
tions (except at i = and 17 = 0, where both methods 
are exact). At finite g/t the two self-energies are differ- 
ent. Moreover, because of the limit d 00, in DMFT 
G itself (not only E) is independent of k. Finally, the 
DMFT evaluation requires self-consistent iterations, and 
is therefore much more involved than that of the MA, 
SCBA and MLF. For these reasons, we do not consider 
the results of the DMFT in the following. 

Interestingly, SCBA also depends on goiui)- Since 

2 

SscsA(t^) ^ Tf^q GscBA{k - q,uj~ ive), we have: 

SsCBa(w) = 5^30 {uJ -UJE - ^SCBA - ^e)) 

= 9^90 {uj-ujE- g^go {uj-2ujE -g^go{uj-3ujE- •■•))) 

On the other hand, for the Holstein model, the MLF 
expression is reminiscent of Eq. ^ ^| : 



G 



MLF 



{k,u;) 



A" 



1 



n! UJ 



Xujf 



n=o -~ - ' — ^' nUJE + lV 

To understand the diagrammatic meaning of the MA 
approximation, we expand Eq. ^ in powers of g^: 

SMA(t^) = g^go{i^ - i^e) + 5'*2go(w - ujE)go{t^ - 2we) 

+9^ [45o(w - wb)5o(^ - 2w£;) -f &gl{uj - uje) 

xgl{uj - 2ujE)goiuJ - 3uje)] + 0{g^) (9) 

showing one contribution of order g^ (which is the correct 
Born expression), 2 of order g*, 10 of order 5^, etc. One 
can verify that this generates the correct total number of 
diagrams in all orders. The difference is that in all MA 
diagrams, each Go(p, f^) free propagator is replaced by a 
momentum averaged 50 (^) function. For example, the 
exact 2'"* order contribution (see Fig. ^: 

^4 
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T72 E ^"(^ - (fi,t^ - ujE)Goik - qi - q2,uj - 2uje) 



[Coik 



qi,uj - uje) + Go{k - q2,uj ~ uje) (10) 
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is replaced within the MA approximation by 

29^(^T.^oiquu;~u;E)j Go(gT>, c.-2c.b)| (H) 

An higher orders are obtained similarly. Let us see why 
this is indeed a good approximation if t ^ g. If i = 0, the 
two expressions are equal. Higher order powers of t come 
from expanding each Go{k, uj) = Go{uj) + egGo(w) + . . . , 
where Go{uj) = {uj + irj)^^. All odd-order powers are 
zero since J2k^^^^ — 0- Consider O(t^) terms in Eq. 

these come either from expanding one of the Go 
to O(t^), in which case they equal their counterparts in 
Eq. 1)11(1 : or they come from 0{t) contributions from two 
different Go lines. In the later case, most terms are zero 
because the two lines generally carry different momenta, 
and J2qi q2 ^91^92 =0. Of 6 such terms generated in Eq. 
((10|l . only one, coming from the outside Go lines of the 
non-crossed diagram, is finite. The error from such terms 
decreases as one goes to higher order diagrams, because 
the percentage of diagrams with one or more pairs of Go 
lines of equal momenta decreases exponentially. Similar 
arguments apply for higher powers in t. It follows that 
MA captures most of the t dependence of each diagram, 
while summing over all diagrams. This suggests that MA 
may be accurate even far from the limit t <^ g. Indeed, 
Eq. clearly shows that MA is also valid for g t. 

For a better idea of the accuracy and range of the MA, 
we consider the sum rules for A{k,uj) = — ilmG(fc, w), 
M„ik) = J^^duJuj''A{k,Lu), which can be evaluated an- 
alytically. The usual approach ^ is based on the equa- 
tions of motion; since MA is also based on them, it should 
fare well. Another approach ^3 is to start with the 
Dyson equation G{k,uj) — Go{k,uj) + Gg(fc, a;)S(A:, -|- 
Gq(A;, uj)T,'^{k, uj) + . . . and the perturbational expansion 
S(fc, uj) = g2s(i) (fc, uj) + .g4s(2) {k,uj) + ..., perform the 
integrals dujuj"G{k, uj) and then take the imaginary 
part. This task is aided by the fact that most terms in 
the integrand decay faster than l/cj as w — *■ oo, and their 
contributions vanish. For n = 0, 1, only Go{k,uj) has fi- 
nite contributions, giving Mo(fc) = 1; Mi{k) — e^. In 
fact. Go contributes an to Mn{k). Next is Ggg^S^^). 
It decays like l/ui^, so it contributes only for n > 2. 
Both SCBA and MA have the exact expression for 
so they both satisfy exactly the n — 2 and 3 sum rules. 
For n > 4, both Ggg^^Cs) and Gg (g^^^^^f contribute. 
Since SCBA ignores one S'^^) diagram, it fails at this 
point, while MA is still exact for n = 4 and 5. MA fails 
at n = 6 because of the approximations in diagrams' ex- 
pressions. Instead of M6(^) = el+g'^[5el + 6t'^{2d'^ -d) + 
4.elujE + 3e|cj| -f Ut^ (e| + 2ekUJE + 2uj%) + 2efew| + uj%] + 
g'^{lMt^ + I2el + 22ekUJE + 25w|) + Ihg^, MA predicts a 
sum rule equal to Me{k) — 2dt^g* (the dimension d enters 
through 2dt^ — ^| ^^'^ higher averages). In other 




g'/(fflgt) g'/(2a)|,t) 



FIG. 2: GS energy, Eos, and quasiparticle weight, Z, for 
uje = 0.5t. d — 1 (left) and d = 2 (right). Inset: line below 
which a second bound peak appears in Id (see text). 

words, MA captures exactly both the dominant power 
in t, e^?, and the dominant power in g, which is or 
ujEg^~^, for even/odd n. This also follows because MA is 
exact both for t = and g = 0, so it can only miss terms 
~ g*t^- These are lost because terms from Gq lines car- 
rying equal momenta are neglected. As discussed, such 
terms are a small minority of all contributions, and in- 
deed MA recovers the vast majority of terms in any M„, 
like in the n = 6 case, showing that it is highly accurate 
not only for t <^ g and t ^ g, but also for intermediary 
values. By contrast, although exact up to n = 3, SCBA 
fails badly at higher n because of the many higher or- 
der diagrams it neglects. For example, for n = 6, SCBA 
predicts 5g^ instead of 15g^ as the leading g term (with 
many 0{g'^) terms missing), showing that SCBA fails for 
g ^ t. Following this analysis we conclude that agree- 
ment with a few sum rules is not meaningful; meaningful 
is to have agreement for the vast majority of terms in all 
sum rules, and in particular for the dominant terms in 
various limits. MA satisfies this restrictive condition. 

MLF also captures both the t = and the g = limits 
exactly. However, this alone does not suffice, either. Di- 
rect evaluation shows that MLF fails badly all sum rules 
with n > 1 if g 7^ 0, as it predicts Mi e~^ek, etc. 

A comparison of MA, SCBA and MLF against GS en- 
ergies and qp weig hts Z = |(GS'|4^J0)|2 obtained with 
DMC [Tl are shown in Fig. ^ for both d = 1 and 2. MA 
compares equally well throughout the Brillouin zone |lCll| . 
The agreement is best for g <C i and g ^ t, but is good 
even for g ^ t. In the inset, we show the line in parame- 
ter space below which a second bound peak appears in Id 
(i.e. the energy of the first excited A: = state satisfies 
El < Eqs +uje)- Above this line, a continuum starts at 
Eas + ^E- The agreement between MA and ED data of 
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FIG. 3: (Color online) A{k,Lj) in Id, for t = l,ujE^ 0.4,7? = 
O.lojE, g = 0.632, 0.894 and 1.265, and fc = and vr. MA data 
(orange line) vs. results from Ref. (black line). 



Ref. is excellent. By contrast, SCBA never predicts a 
second bound state (it always finds the continuum) ; MLF 
always predicts a peak at Egs +^e, never a continuum. 

A comparison of the MA spectral weight A{k,uj) (or- 
ange line) in Id, for k — and n/a and three values of 
g is shown in Fig. 01 against data (black line) based on 
a variational method [Sj. The agreement is again decent, 
especially since g ~ t in all three cases. Other compar- 
isons with numerical data are of similar quality jlCll |. 

MA thus provides an accurate, yet simple and fast way 
to study G{k,u}) for Holstein polarons in any dimension 
for any e^. The question is whether this approach can be 
expanded beyond the Holstein Hamiltonian. While this 
issue is still under investigation jlQj . one possible route 
is provided by this generalization of Eq. : 



t{k,uj) 



-y- 



Ei{k,qi,\i) 
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E2{k, (fi, 92, Al, A2) 



(12) 
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Eiik, qi,\i) = l5fc^g-iAi PG'o(fci, t^i), E2{k, qi,q2, Xi, X2) 



where ki ^ k — ql, ki 2 = k — qi — q2 and Ui — uj — uj^.^x., 
^1,2 — uj — ^qi,Xi — ^q.2.\-2- This describes coupling to 
several branches of bosons with dispersions ijJ^\^ with a 
vertex if k is the momentum of the incoming parti- 
cle and qX is the momentum and branch of the emitted 
boson. While not exact, S(A:,ti;) is much more accurate 
than Ea/^i: it gives aU V\ 2"'^ and 9 out of 10 of the S'''* 
order diagrams exactly. The wrong S' '^ order diagram has 
one of its five Go lines with a wrong momentum and en- 
ergy. All higher order diagrams' numbers and topologies 
are generated correctly, a small fraction of them having 
some mislabelled Gq lines. One can trace the first failing 
of a sum rule, due to the wrong 3'"'' order diagram, to now 



occur in M^ik). For simplicity, let us assume that boson 
frequencies uj\ and vertices g\ depend only on the branch. 
Then, in Mg,{k), a term (|s'aP)(wa|3aP)^ is replaced by 



, and an extra term 2dt^ {\g xi^)^ is gener 



{\9>.??{^l\9>.?) 
ated, but all other terms including the dominant g terms 
lOSd^Ap)"* are exact (here, (/) — /). Using a further 
MA approximation removes the need to evaluate the mo- 
mentum sums in (|12f) by replacing all Gq{p, fi) go{^)] 
this results in an error in MQ{k) (a missing 2dt'^{\g\\'^)'^) 
but would speed up calculations significantly and, as for 
the Holstein model, should have a limited effect on ac- 
curacy. Further simplifications are possible if the boson 
frequencies are close to each other. One can show that 
if i = and aU lux = uJe, Eq. Q with g"^ = {\gx\^) is 
the exact self-energy [lO| . We do not know if this iden- 
tity has been noted before. For close-by phonon energies, 
one can then also remove the branch sums in H12|) by us- 
ing Eq. C3 with g2 ^ (|5Ap), i^E = (wA|ffAp)/5^ in 
which case one gets an error starting with Mi{k), where 
(w^l^Ap) — >■ 9^^%, but all dominant terms are correct 
in all orders. Eq. (|12|l may also be easy to estimate for 
highly anisotropic g^-, in which case the BZ sums reduce 
to summations over a few hot spots. Of course, tests 
against numerical results are needed to verify all this. 

To conclude, progress has been made in a very old 
problem, by finding a simple yet highly accurate approx- 
imation for G{k,Lu) of the Holstein polaron. A path to 
possibly more exciting results has also been uncovered. 
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